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Abstract 

A CFD method for solving combustor flowfields at all 
speeds on complex configurations is presented. The ap- 
proach is based on the ALLSPD-3D code which uses the 
compressible formulation of the flow equations including 
real gas effects, nonequilibrium chemistry and spray com- 
bustion. To facilitate the analysis of complex geometries, 
the chimera grid method is utilized. To the best of our 
knowledge, this is the first application of the chimera 
scheme to reacting flows. In order to evaluate the effec- 
tiveness of this numerical approach, several benchmark 
calculations of subsonic flows are presented. These in- 
clude steady and unsteady flows, and bluff-body stabi- 
lized spray and premixed combustion flames. The results 
demonstrate the effectiveness of the combined ALLSPD- 
3D/chimera grid approach for analyzing subsonic com- 
bustor flowfields down to the incompressible limit. 

1. Introduction 

The numerical simulation of flow within a com- 
bustor requires a CFD code capable of solving the 
three-dimensional Navier-Stokes equations with finite- 
rate chemistry on complex geometries. The compressible 
formulation is usually required since in many combustor 
applications velocities can range from low subsonic to su- 
personic. Even when the Mach number is low throughout 
the combustor, a compressible formulation is still neces- 
sary because of the wide density variation resulting from 
the combustion process. The computation of flowfields 
involving regions of low subsonic velocities using a com- 
pressible formulation is complicated by two main numeri- 
cal difficulties. First, the flow equations become “stiff” at 
low Mach numbers. The degree of stiffness is determined 
by the ratio of the largest to the smallest eigenvalues. 
When the Mach number approaches zero, this ratio be- 
comes infinitely large, and the Navier-Stokes equations 
become singular. In theory this difficulty can be circum- 
vented by using implicit schemes. However, due to fac- 
torization errors, such schemes usually produce very poor 
convergence rates for problems involving very low Mach 
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numbers. Second, the pressure gradient terms in the mo- 
mentum equations become singular as the Mach number 
approaches zero. (In nondimensional variables, the pres- 
sure term is of order 1/M 2 while the convective term is 
of order unity). This singularity problem causes large 
roundoff errors at low Mach numbers. 

At NASA Lewis Research Center a new CFD code, 
ALLSPD-3D 1 , is being developed for efficiently solving 
combustor flowfields at all speeds. This code solves the 
three-dimensional chemical nonequilibrium Navier-Stokes 
equations over a wide range of Mach numbers. The 
method used is based on the strong conservation form 
of the governing equations, but utilizes primitive vari- 
ables as unknowns. To overcome the difficulties for low 
Mach number calculations a preconditioning matrix is 
introduced 2 . The singular behavior of the pressure gra- 
dient terms is circumvented by expressing the pressure as 
a sum of a reference pressure and a gauge pressure 2 , as 
discussed below. Real gas properties and nonequilibrium 
chemistry are considered. Turbulence is modeled with a 
low-Reynolds number k - e model developed by Shih et 
al. 3 , and a stochastic model is used to compute liquid 
spray combustion 4 . 

In addition to having efficient numerics and accurate 
physical models, the computational method must be able 
to handle complex geometries. An efficient approach for 
this purpose is the domain decomposition method known 
as the chimera scheme 5 . This approach uses a compos- 
ite of independent overlapping grids to discretize a com- 
plex geometrical configuration. The main advantages of 
the chimera method include; 1) simplification of the grid 
generation process by permitting a modular approach; 
2) improvement in computational efficiency by provid- 
ing better flowfield resolution with fewer grid points; 3) 
simplification in the application of boundary conditions; 
and 4) the possibility of allowing relative motion between 
system components without requiring a new grid to be 
generated at each time step. 

Some of the disadvantages to the use of chimera grids 
are the complexity in the bookkeeping that tracks rela- 
tionships among grids, and the fact that conservation 
is not strictly enforced at interface boundaries. In the 
present work we show that for low Mach number flows, 
mass conservation does not appear to present a serious 
problem. 

The objective of this study is to evaluate the effective- 
ness of the combined ALLSPD-3D/chimera grid scheme 
approach for computing combustor flowfields. To the best 
of our knowledge, the chimera grid scheme has not been 
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applied to reacting flows previously. The numerical for- 
mulation is described below, followed by a description of 
the numerical method and chimera scheme. Results are 
then presented for several benchmark tests including a 
channel flow and steady and unsteady flow over a bluff- 
body. Results are also presented for a bluff-body stabi- 
lized spray combustion flame. 

2. Mathematical Formulation 

The three-dimensional, unsteady, compressible, 
density-weighted time-averaged Navier-Stokes equations 
and species transport equations for a chemically react- 
ing gas of N species written in generalized nonorthogonal 
coordinates can be expressed as 



= Hc + H/, (1) 

where the vectors Q, E, F,G, E v , F Vl G V} H c and H/ 
are defined as 



E = j(**Q + €,E + *,F + &G), 

E = j (*7tQ + + rjyF + i/ z G), 

G = j(<*Q + <.E + <,F + GG), 

E» = y(£*E v -h^yFi; + ^Gt,), 

E v = —(rj r E u + *}yF v 4* ^Gu), 

G v = ~ (CrEt, *f £ y F v 4- 0G t) ), 

H c = jH c , 

H, = I*. 

In the above expressions, r, £, r) and ( are the temporal 
and spatial generalized coordinates and £ t , T) t and are 
the grid speed terms. r)x, Vy, Vz, Cr, and 

C* are the metric terms and J is the transformation Ja- 
cobian. The vectors Q, E, F, G, E„, F„ and G„ in the 
above definitions are 

Q = (p,pu,pv,pw,pE,pK,pt,pYi,- ■ ■,pYs-y) T , 

E = (p«,pu 2 +p,puv,puw,(pE + p)u,puK,puc, 

puY u ■ ■ ■,puY N _ 1 ) T , 

F = (pv,puv,pv 7 +p, pvw, (pE + p)v, pVK, pvc, 
pvYi , ■ ■ -,pvY N _ i) T , 

G = (pw,puw,pvw,pw 2 + p,(pE + p)w,pwK, 
pwt,pwYi, ■ ■ -,pwY N - ,) T . 

E t) = (0, T XXl T X y, T xt , UT XX + VT X y + WT X . + q X(t 

\T 

r XK I Txt } 9rl , * ■ 0x;v-l ) J 


E V (0, Ty X , Tyy , Ty Z , UTy j; + VTyy + XV Ty Z + Qy ^ , 

r y*i r y«i tfyi > ' • ' J , 

Gf = (0, r lx ,T Z y,T iz ,uT zx + VT Z y + wt. l + q ltt 

T ZKi T ZC, Qil, * * ' ) , 

where p } p, u, v, iu, «, c and Y, represent the density, pres- 
sure, Cartesian velocity components, turbulent kinetic 
energy, dissipation rate of turbulent kinetic energy and 
species mass fraction , respectively; E — e+^(u 2 -fv 2 -btf 2 ) 
is the total internal energy with e being the thermody- 
namic internal energy; The normal and shear stresses 
( r n, r xy ,...), energy (g re ,...), species and tur- 

bulent diffusion (r x *,...) fluxes are given in Chen and 
Shuen 4 . 

The first term in Eq. (1), r<9Q/<9r*, is the time pre- 
conditioning term which is added to overcome the difficul- 
ties encountered in low Mach number calculations when 
the compressible formulation is used. The detailed theory 
about this specific treatment has been described in Shuen 
et al. 2 . The expressions in this term are 
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where Q is the primitive dependent variable (to be 
solved), r is the preconditioning matrix, r m is the pseudo- 
time and h — e - h ^ is the specific enthalpy of the gas 
mixture. /? is a scaling factor and is taken to be 

/? = u 2 + v 2 + w 7 . 

As described previously, this preconditioning technique 
allows us to calculate very low Mach number flows with- 
out difficulties. To circumvent the pressure gradient sin- 
gularity problem, the pressure has been expressed as a 
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sum of a constant reference pressure part, p 0 , and a gauge 
pressure part, p g : 

P = P0+Pg 

The reference pressure is usually selected as the 
freestream pressure. When this expression for the pres- 
sure is used in Eq. (1), the pressure gradient term in the 
nondimensional momentum equations becomes of order 
one as the Mach number approaches zero, and the singu- 
larity is eliminated from the equations. 

The source term vectors H c and H / are 


In reacting flow calculations, the evaluation of thermo- 
physical properties is of vital importance. In this paper, 
the values of Cp,, fc/, , and pn for each species are deter- 
mined by fourth-order polynomials of temperature 7 . The 
specific heat of the gas mixture is obtained by mass con- 
centration weighting of individual species. The thermal 
conductivity and viscosity of the mixture, however, are 
calculated using Wilke’s mixing rule 8 . The binary mass 
diffusivity Dij between species i and j is obtained using 
the Chapman-Enskog theory 8 . 
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The quantities related to the source term in the tur- 
bulent equations (®, ci, etc.) are given by Chen and 
Shuen 4 . The term 5, is the chemical source term for 
species i, and is evaluated using standard methods (see 
for example Ref. 6). The spray source term H /, spray 
model, liquid-phase Lagrangian equations of motion and 
droplet mass and heat transfer equations are described in 
detail by Chen and Shuen 4 and will not be repeated here. 
The temperature and density are calculated iteratively 
from the following equations 


3. Numerical Method 
3.1 Iterative Scheme 

The gas-phase governing equations (Eq. (1)), which in- 
clude the fluid dynamics, turbulence, and species conser- 
vation equations, could in principle be solved in a fully 
coupled manner. This approach would yield the best sta- 
bility, robustness and convergence rate. However, such an 
approach requires large amounts of memory and expen- 
sive matrix inversions, which translates into more CPU 
time per iteration. This is specially true when many 
species are considered in the analysis. In order to com- 
promise between the degree of coupling and an efficient 
use of computer resources a partially decoupled approach 
was adopted for solving Eq. (1). The governing equations 
are divided into three sets, i.e., five flow equations, two 
K-e turbulent equations and N — 1 species equations as 
discussed in Ref. 1. Each individual set is solved in a 
coupled manner and iteratively. 

The numerical scheme used in this study, and applied 
to each of the three sets of equations, is the LUSGS fully 
implicit method which is described in detail by Chen and 
Shuen 1,4 and, for high-speed applications, by Yungster 
and Radhakrishnan 9 . The discretized form of Eq. (1) is 
expressed as follows. 

{i 4 {r-Ar*[S + T + (^-£R«£) 

+(^T ~ £; Rr w£;) + (t^ - ^ R cc^)] 

+D 2 }P + (1 - u)I}AQ = { — Ar’(R) p }ij, (2) 

where 

( R )p = + idijj-) + 

gL G-G.) _ _d c + D „ (3) 


hi = h> u + j T C Pi dT, 

JT re j 


N 

p — pRuT ^ 


1=1 


Wi' 


where Ru and T re f are the universal gas constant and 
reference temperature for thermodynamic properties, and 
Wi , C Pt) h { , h J. are the molecular weight, constant pres- 
sure specific heat, thermodynamic enthalpy, and heat of 
formation of species i, respectively. 


where the superscript p denotes the previous iteration 
level, S and T are the Jacobian matrices for chemical 
and turbulent source terms, respectively, A, B and C 
are the inviscid term Jacobians and R^, and R^ 
are the viscous term Jacobians. The expressions for these 
Jacobians are given in Chen et al 1 . Note that the physical 
time term has not been included in the above expression, 
and therefore the solution marches forward in pseudo- 
time. The second and forth order dissipation terms for 
LUSGS scheme are defined as 

D 2 = -0.5a [^7 ( |A>t I r) + ^7(|A B |r) + |A c |r)], 
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D 4 - o-r( |A>i| + |Ab| + | Ac | 

where lA^I, |Ab| and |Ac| are the absolute value of the 
maximum eigenvalue in the £ , rj and £ directions, respec- 
tively. The D c in Eq. (3) represents a source term due 
to the partially decoupled procedure 1 . The a and a are 
defined in Chen and Shuen 4 . The term it is a flag used 
with the chimera method which is described below. 

3.2 Chimera method 

In the chimera scheme approach, a system of relatively 
simple overlapping grids, each describing a component of 
a complex configuration, is combined into a composite 
grid to yield solutions for complex flow fields. A method 
to interconnect the grids is needed to determine when 
points of one grid fall within a body boundary of an- 
other (grid hole points), and to supply pointers so that 
one grid can provide data to update the boundary points 
of another. In this study, the code PEGSUS 10 has been 
used for this purpose. The information is passed from one 
grid to another via trilinear interpolation. In order to dis- 
tinguish the hole and grid boundary points from regular 
computational points, a blanking array, it , is used in the 
flow solver. For a hole or grid boundary point, it is set 
to zero; otherwise it is set to one. In order to exclude the 
hole and grid boundary points from the solution proce- 
dure, all the terms in the right hand side and left hand 
side of equation (2) are multiplied by it , and the value of 
(1 — it ) I is added to the left hand side of this equation. 
Thus, tor a regular computational point for which it — 1, 
the normal scheme is maintained, but when it = 0 the 
scheme reduces to AQ = 0 and thus Q is not changed at 
a hole or boundary point. The solutions, Q, correspond- 
ing to the boundary points with it = 0 are updated by 
the trilinear interpolation as described above. 

4. Results 

Results using the ALLSPD-3D/chimera approach are 
presented for four subsonic benchmark test cases: 1) non- 
reacting laminar channel flow; 2) nonreacting steady and 
unsteady bluff-body flows; 3) bluff-body stabilized spray 
combustion flame; and 4) bluff-body stabilized premixed 
combustion flame. 

4.1 Laminar Channel Flow 

The first case considers a low speed (V=0.051 m/sec), 
laminar flow in a channel with an arc bump. In oraer to 
test the numerical method, solutions are computed on a 
single zone grid (41x21x17) and a 2-block chimera grid 
with block dimensions of (27x21x171 and (20x23x19). 
Figure 1 shows the single zone and the 2-block chimera 
grids. For the chimera grid, the overlapped region be- 
tween two grids is indicated in Fig. lb. The boundary 
conditions are (see Fig. la): £min '• subsonic inflow; 
subsonic outflow; 7? min * no-slip wall; T] max : symmetry; 
Cnin : symmetry; Cmax- no-slip wall. 

The convergence history is presented in Fig. 2. It is 
noticed that the rate of convergence for the chimera grid 
is nearly a s good as that for the single zone grid. A four 
order of magnitude reduction in the residual is obtained 


in approximately 1000 iterations. It is also worth pointing 
out that the overhead in memory requirements and CPU 
time due to the use of chimera grids is negligible. 

The velocity magnitude along the center of the channel 
is plotted in Fig. 3. Excellent agreement between the two 
solutions is obtained. 

4.2 Laminar Steady and Unsteady Bluff-Body 
Flows 

In this case, a circular cylinder is placed inside the 
channel described in the previous section. Individual 
overlapping grids are generated for the channel and cylin- 
der as shown in Fig. 4. Laminar results were obtained 
for two inflow velocities corresponding to Reynolds num- 
bers (based on the cylinder diameter) of Rej= 36.7 and 
Re<i=73. 4. It is known 11 that for an unbound flow over a 
circular cylinder, the upper limit of Reynolds number for 
a steady-state solution is about 40. In the present case, 
the flowfield around the cylinder is constrained by the 
channel walls, however, this is not expected to change dra- 
matically the limiting Reynolds number for steady flow. 
Indeed, the results indicate that for the low Reynolds 
number case (/?e < j= 36.7) the flow reaches a steady state, 
while for the /Jej=73.4 case the flow exhibits an unsteady, 
periodic vortex shedding behavior. 

The steady-state solution for the Rtd= 36.7 case is pre- 
sented in Fig. 5 in the form of velocity vectors. One of the 
main concerns regarding the use of the chimera method 
is the fact that since interpolation is used to connect 
grids, conservation is not strictly enforced at the inter- 
face boundaries. To test the overall mass conservation, 
the mass flow rate across the various £ = const channel 
planes was computed for this case, and the results plotted 
in Fig. 6. The gap in the plot corresponds to the loca- 
tion occupied by the circular cylinder. This plot clearly 
indicates that the error due to interpolation is very small, 
less than 0.5%. 

Figure 7 shows a series of velocity vectors at selected 
pseudo-times during one of the periodic cycles of the 
J Rc (i =73.4 case. This figure illustrates the changing pat- 
terns of the wake. In Fig. 7f, the vortex structure behind 
the cylinder is nearly identical to that observed at the 
beginning of the cycle (Fig. 7a). The unsteady periodic 
nature of this flowfield can also be clearly observed in the 
convergence plot discussed later. 

4.3 Nonpremixed Bluff-Body Stabilized Flame 

The application of the present chimera grid method to 
combustion flow is demonstrated here using the same grid 
setup as the previous case. The circular cylinder is used 
as a flame holding device to stabilize the flame. For the 
present test purpose, the flow is assumed laminar and the 
Reynolds number based on the cylinder diameter and in- 
let air velocity is 57.7. The inlet air temperature is 600 
K. Liquid methanol was injected from 10 injection holes 
behind the cylinder. The schematic of this injection is 
depicted in Fig. 8. The liquid methanol was assumed to 
be fully atomized with the initial drop diameters ranging 
from 20 to 100 fi m. Five species (CH 3 OH , O 2 , ^2i CO 2 , 
and H 2 O) were considered in this calculation. A single- 
step global reaction chemistry model reported in West- 
brook and Dryer 12 was used. The fuel/air ratio is 0.1. 
The results of this calculation at three selected stream- 
wise planes are presented in Fig. 9. Contour plots for fuel 
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and oxidizer mass fractions are shown in Figs. 9a and 9b, 
and temperature contours are shown in Fig. 9c. Near 
the liquid injection, there is a relatively fuel-rich core re- 
gion due to significant amount of evaporation of the liquid 
methanol. This can be clearly seen in the fuel mass frac- 
tion contours. The flame was established along the outer 
edge of this fuel-rich core region where gaseous fuel and 
oxidizer are well mixed. The first temperature contour 
plane behind the cylinder indicates this flame pattern. 
The second contour plane indicates that at this station 
all the fuel has been consumed, and a nearly circular high 
temperature region is established. Further downstream, 
the surrounding air mixes with the combustion products 
reducing the temperature levels significantly. 

For the previous three bluff-body non-reacting and re- 
acting flows, the convergence histories are shown in Fig. 
10. For non-reacting flows, steady state convergence is ob- 
served for the low Reynolds case {Rtd— 36.7), while a peri- 
odic convergence pattern exhibits for the higher Reynolds 
number flow (Re<t=73A). This periodic convergence pat- 
tern is consistent with the vortex shedding flow previously 
discussed. For the combustion convergence curve, a spike 
near 100 iterations indicates the time of ignition. Several 
small peaks near the end of the convergence curve are due 
to the numerical treatment of the spray source terms that 
are updated every 20 iterations in the present calculation. 

4.4 Premixed Bluff-Body Stabilized Flame 

To validate the current ALLSPD-3D/chimera method 
on the accuracy of the combustion flow predictions, a pre- 
mixed propane flame with available experimental and nu- 
merical data 13 was chosen in this study. The test case 
consists of a straight channel with a rectangular cross 
section and a wedge flame holder. Although the exper- 
imental test was conducted in a fully three-dimensional 
rig, a two-dimensional flow is calculated along the sym- 
metry plane where the experimental data and other nu- 
merical predictions are available to be compared with. A 
schematic, taken from Ref. 13, of the experimental test 
setup is shown in Fig. 11. The gaseous propane fuel is 
premixed with air before entering the combustion sec- 
tion. The inlet temperature and velocity are 288 K and 
17.06 m/s, respectively. The equivalence ratio is 0.61. 
Five species (CaJ/g, 0 2 , C0 2 , and H 2 0) were con- 

sidered in this calculation. A single-step global reaction 
chemistry model reported in Westbrook and Dryer 12 was 
used. 

Two grids were used in this study. Figure 12 shows a 
Cartesian grid (175x41) for the channel and an O-type 
grid (31x130) around the wedge flame holder. Since the 
combustor is symmetric in the y-direction, only half of 
the domain was actually calculated. The converged re- 
sults are shown in Fig. 13. The temperature contours are 
presented in Fig. 13a. The combustion flame is stablized 
behind the flame holder with a vivid flame front that 
starts near the tip of the flame holder and extends further 
downstream. Figures 13b and 13c show the comparison 
for temperature profiles between the current predictions 
and two other data at two streamwise stations behind the 
flame holder. The first station is located 0.15 meter be- 
hind the flame holder and the second station is at 0.35 
meter. The temperature profile at the 0.15 meter station 
compares favorably with the numerical predictions and 
experimental data reported in Ref. 13. However, both 


current results and the other numerical predictions at the 

0.35 meter station do not agree that well with the experi- 
mental data. As indicated in Ref. 13, the diffusion /mixing 
rate is too slow compared to the experimental data. This 
can be seen in both temperature profiles. The isotropic 
assumption of the current low Reynolds turbulence model 
may not be adequate in predicting the non-isotropic tur- 
bulence in the shear layer region close to the flame holder. 
Although many issues exist and further study is needed to 
improve the accuracy of the current simulation, the use of 
the chimera method presented no difficulty in obtaining 
reasonably good results 

5. Conclusions 

The successful implementation of the chimera grid 
scheme into the ALLSPD-3D code has extended its ca- 
pability to handle complex configurations. The use of 
the chimera method for domain decomposition resulted 
in very low overhead in CPU time and memory require- 
ments, and did not degrade the excellent convergence 
characteristics of the numerical method. The effective- 
ness of this numerical approach was demonstrated by 
computing several subsonic, nonreacting, steady and un- 
steady flows, and bluff-body stabilized spray and pre- 
mixed combustion flames. Mass conservation error was 
shown to be very small (< 0.5%) for these cases. Also, 
the large gradients in the flow variables (i.e., temperature, 
density, species concentrations, etc) resulting from com- 
bustion processes did not have a negative impact on the 
interpolation process. The results of this work indicate 
that the combined ALLSPD-3D/chimera grid approach is 
a powerful tool for computing subsonic combustor flow- 
fields down to the incompressible limit. 
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Figure 2. Convergence history 



channel flow, (a) single zone grid; (b) 2-block x (m) 

chimera grid. 

Figure 3. Comparison of velocity magnitude 
for single zone and chimera grids along chan- 
nel centerline. 
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Figure 4. Two-zone chimera grid for channel/cylinder problem. 
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Figure 8. Schematic of liquid methanol spray injection. 
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Figure 9. Spray combustion Howlield; (a) methanol mass fraction; (b) oxygen mass fraction: 
tc) temperature (K). 
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Figure 11. Schematic of premixed combustion experiment with bluff body flameholder (Ref. 13). 



Figure 12. Two-zone chimera grid. 
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Figure 13. (a) Temperature contours and (b) temperature distribution along the x=0.15 m and 
x=0.35 m stations for the premixed combustor. 
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